home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / VECTOR / MATH2.PAS < prev   
Encoding:
Pascal/Delphi Source File  |  1994-07-07  |  12.9 KB  |  655 lines
  1. unit Math2;
  2.  
  3. interface
  4.  
  5. uses winprocs,math1;
  6.  
  7.  
  8. FUNCTION EVEN(TEST,YES,NO:REAL):REAL;
  9. FUNCTION POWER1(Y,X,DEFAULT:REAL):REAL;
  10. FUNCTION QUAD(Y,X:REAL):INTEGER;
  11. FUNCTION ARCSIN(S:REAL):REAL;
  12. FUNCTION ARCTAN1(NUM,DENOM:REAL):REAL;
  13. FUNCTION VLENGTH(V1,V2,V3,E1,E2,E3:REAL):REAL;
  14. FUNCTION LOC(B,C,A:REAL;MODE:INTEGER):REAL;
  15. FUNCTION LOS(A,AANGLE,B,DEFAULT:REAL;MODE:INTEGER):REAL;
  16. PROCEDURE COMPLEXCONVERT(A1,A2:REAL;MODE:INTEGER;VAR O);
  17. PROCEDURE SPHERICAL_XYZ(S1,S2,S3:REAL;MODE:INTEGER;VAR SA);
  18. PROCEDURE XYZROTTRANS(X,Y,Z,X1,Y1,Z1,X2,Y2,Z2,T1,T2,T3,XM,YM,ZM:REAL;MODE:INTEGER;VAR R);
  19. PROCEDURE BAIRSTOW(VAR A,B;EPSILON:REAL;IT1,DEGREE:INTEGER);
  20. PROCEDURE MULTIPLYCR(X,Y,X1,Y1:REAL;VAR Z);
  21. PROCEDURE DEVIDECR(N1,N2,D1,D2:REAL;VAR R);
  22. PROCEDURE ADDSUBCP(A1,A2,B1,B2:REAL;MODE:INTEGER;VAR R);
  23. PROCEDURE POWERCR_CP(A1,A2,B1,B2:REAL;MODE:INTEGER;VAR R);
  24. PROCEDURE RAYS_XYZ(X,Y,Z,X1,Y1,Z1,X2,Y2,Z2,XM1,YM1,ZM1:REAL;COUNT:INTEGER;VAR RIN,ROUT);
  25.  
  26. implementation
  27.  
  28. FUNCTION EVEN(TEST,YES,NO:REAL):REAL;
  29. BEGIN
  30. IF TEST/2=TRUNC(TEST/2) THEN
  31. BEGIN
  32. EVEN:=YES;
  33. END
  34. ELSE
  35. BEGIN
  36. EVEN:=NO;
  37. END;
  38. END;
  39.  
  40. FUNCTION POWER1(Y,X,DEFAULT:REAL):REAL;
  41. VAR
  42. Y1:INTEGER;
  43. BEGIN
  44. IF Y<=0 THEN
  45. BEGIN
  46. Y1:=-1;
  47. IF Y=0 THEN
  48. BEGIN
  49. POWER1:=0;
  50. END
  51. ELSE
  52. BEGIN
  53. IF TRUNC(X)=X THEN
  54. BEGIN
  55. POWER1:=(Y1+EVEN(X,2,0))*(EXP(X*LN(ABS(Y))));
  56. END
  57. ELSE
  58. BEGIN
  59. POWER1:=DEFAULT;
  60. END;
  61. END;
  62. END
  63. ELSE
  64. BEGIN
  65. POWER1:=EXP(X*LN(Y));
  66. END;
  67. END;
  68.  
  69. FUNCTION QUAD(Y,X:REAL):INTEGER;
  70. BEGIN
  71. {quad is used with arctan1.  assigns values from 0 to 8 as you rotate in x-y plane
  72. around z axis.  numbers represent axis or quadrants}
  73. IF (Y=0)AND(X=0) THEN
  74. BEGIN
  75. QUAD:=0;
  76. END
  77. ELSE
  78. BEGIN
  79. IF (Y=0)OR(X=0) THEN
  80. BEGIN
  81. IF Y=0 THEN
  82. BEGIN
  83. IF X<0 THEN
  84. BEGIN
  85. QUAD:=5;
  86. END
  87. ELSE
  88. BEGIN
  89. QUAD:=1;
  90. END;
  91. END
  92. ELSE
  93. BEGIN
  94. IF Y<0 THEN
  95. BEGIN
  96. QUAD:=7;
  97. END
  98. ELSE
  99. BEGIN
  100. QUAD:=3;
  101. END;
  102. END;
  103. END
  104. ELSE
  105. BEGIN
  106. IF Y>0 THEN
  107. BEGIN
  108. IF X>0 THEN
  109. BEGIN
  110. QUAD:=2;
  111. END
  112. ELSE
  113. BEGIN
  114. QUAD:=4;
  115. END;
  116. END
  117. ELSE
  118. BEGIN
  119. IF X>0 THEN
  120. BEGIN
  121. QUAD:=8;
  122. END
  123. ELSE
  124. BEGIN
  125. QUAD:=6;
  126. END;
  127. END;
  128. END;
  129. END;
  130. END;
  131.  
  132. FUNCTION ARCSIN(S:REAL):REAL;
  133. VAR
  134. C,F:REAL;
  135. BEGIN
  136. IF ABS(S)>1 THEN
  137. BEGIN
  138. ARCSIN:=0;
  139. END
  140. ELSE
  141. BEGIN
  142. IF ABS(S)=1 THEN
  143. BEGIN
  144. ARCSIN:=(S/ABS(S))*PI/2;
  145. END
  146. ELSE
  147. BEGIN
  148. C:=SQRT(1-SQR(S));
  149. F:=ARCTAN(S/ABS(C));
  150. ARCSIN:=F;
  151. END;
  152. END;
  153. END;
  154.  
  155. FUNCTION ARCTAN1(NUM,DENOM:REAL):REAL;
  156. VAR
  157. Q:INTEGER;
  158. BEGIN
  159. {with num as numerator and denom as denominator, arctan1 gives values from 0 to 2*pi
  160. instead of -pi/2 to pi/2, good for use with rotating object}
  161. Q:=QUAD(NUM,DENOM);
  162. CASE Q OF
  163. 0:ARCTAN1:=0;
  164. 1:ARCTAN1:=0;
  165. 2:ARCTAN1:=ARCTAN(NUM/DENOM);
  166. 3:ARCTAN1:=PI/2;
  167. 4:ARCTAN1:=ARCTAN(NUM/DENOM)+PI;
  168. 5:ARCTAN1:=PI;
  169. 6:ARCTAN1:=ARCTAN(NUM/DENOM)+PI;
  170. 7:ARCTAN1:=3*PI/2;
  171. 8:ARCTAN1:=ARCTAN(NUM/DENOM)+2*PI;
  172. ELSE
  173. ARCTAN1:=0;
  174. END;
  175. END;
  176.  
  177. FUNCTION VLENGTH(V1,V2,V3,E1,E2,E3:REAL):REAL;
  178. BEGIN
  179. {gives vector length for head of vector v1v2v3 to tail e1e2e3}
  180. VLENGTH:=SQRT(((v1-e1)*(V1-E1))+((v2-e2)*(V2-E2))+((v3-e3)*(V3-E3)));
  181. END;
  182.  
  183. FUNCTION LOC(B,C,A:REAL;MODE:INTEGER):REAL;
  184. BEGIN
  185. {law if cosines, b and c are sides, a can be angle or side, for mode=1, a angle, result is length
  186. of third side opposite a,  for mode<>1 then a,b,c, sides and angle opposite a is solved for}
  187. IF MODE=1 THEN
  188. BEGIN
  189. IF (B=0)OR(C=0) THEN
  190. BEGIN
  191. LOC:=SQRT(SQR(B)+SQR(C));
  192. END
  193. ELSE
  194. BEGIN
  195. LOC:=SQRT(SQR(B)+SQR(C)-(2*B*C*COS(A)));
  196. END;
  197. END
  198. ELSE
  199. BEGIN
  200. IF (B=0)OR(C=0) THEN
  201. BEGIN
  202. LOC:=0;
  203. END
  204. ELSE
  205. BEGIN
  206. LOC:=ARCCOS((SQR(B)+SQR(C)-SQR(A))/(2*B*C));
  207. END;
  208. END;
  209. END;
  210.  
  211. FUNCTION LOS(A,AANGLE,B,DEFAULT:REAL;MODE:INTEGER):REAL;
  212. BEGIN
  213. {law of sines, aangle opposite a side, mode 1,a side, b angle, solves for side opposite b angle
  214. mode<>1, a,b sides solves for angle opposite to b} 
  215. IF MODE=1 THEN
  216. BEGIN
  217. IF (A=0)OR(AANGLE=0) THEN
  218. BEGIN
  219. LOS:=DEFAULT;
  220. END
  221. ELSE
  222. BEGIN
  223. LOS:=(A*SIN(B))/(SIN(AANGLE));
  224. END;
  225. END
  226. ELSE
  227. BEGIN
  228. IF (A=0)OR(AANGLE=0) THEN
  229. BEGIN
  230. LOS:=DEFAULT;
  231. END
  232. ELSE
  233. BEGIN
  234. LOS:=ARCSIN((B*SIN(AANGLE))/A);
  235. END;
  236. END;
  237. END;
  238.  
  239. PROCEDURE COMPLEXCONVERT(A1,A2:REAL;MODE:INTEGER;VAR O);
  240. TYPE
  241. O1=ARRAY[1..2] OF REAL;
  242. BEGIN
  243. {mode 1:a1 magnitude a2 angle, converts to rectangular coordinates and places them in o:array[1..2] of real
  244. mode<>1:a1 real a2 imaginary converts to polar and places magnitude in o[1] and angle in o[2]}
  245. IF MODE=1 THEN
  246. BEGIN
  247. O1(O)[1]:=A1*COS(A2);
  248. O1(O)[2]:=A1*SIN(A2);
  249. END
  250. ELSE
  251. BEGIN
  252. O1(O)[1]:=SQRT(SQR(A1)+SQR(A2));
  253. O1(O)[2]:=ARCTAN1(A2,A1);
  254. END;
  255. END;
  256.  
  257.  
  258. PROCEDURE SPHERICAL_XYZ(S1,S2,S3:REAL;MODE:INTEGER;VAR SA);
  259. TYPE
  260. SXYZ=ARRAY[1..3] OF REAL;
  261. VAR
  262. S:REAL;
  263. BEGIN
  264. {in mode 1:s1=ro,s2=theta,s3=phi converts to xyz and places them in sa:array[1..3] of real
  265. {in mode<>1 s1=x,s2=y,s3:=z converts to ro,theta,phi and places them in sa in that order}
  266. IF MODE=1 THEN
  267. BEGIN
  268. SXYZ(SA)[1]:=S1*COS(S2)*SIN(S3);
  269. SXYZ(SA)[2]:=S1*SIN(S2)*SIN(S3);
  270. SXYZ(SA)[3]:=S1*COS(S3);
  271. END
  272. ELSE
  273. BEGIN
  274. S:=SQRT(SQR(S1)+SQR(S2)+SQR(S3));
  275. SXYZ(SA)[1]:=S;
  276. SXYZ(SA)[2]:=ARCTAN1(S2,S1);
  277. SXYZ(SA)[3]:=ARCCOS(S3/S);
  278. END;
  279. END;
  280.  
  281. PROCEDURE XYZROTTRANS(X,Y,Z,X1,Y1,Z1,X2,Y2,Z2,T1,T2,T3,XM,YM,ZM:REAL;MODE:INTEGER;VAR R);
  282. LABEL 3,4,5;
  283. TYPE
  284. R1=ARRAY[1..3] OF REAL;
  285. VAR
  286. U:ARRAY[1..3] OF REAL;
  287. U1:ARRAY[1..3] OF REAL;
  288. U2:ARRAY[1..3] OF REAL;
  289. U3:ARRAY[1..26] OF REAL;
  290. U4:ARRAY[1..10] OF REAL;
  291. H10:ARRAY[1..3,1..4] OF REAL;
  292. M,AXIS:INTEGER;
  293. A,B,LEN,XX,XY,XZ,ZX,ZY,ZZ,YX,YY,YZ,XN,XN1,XN2,XN3,YN,YN1,YN2,YN3,
  294. ZN,ZN1,ZN2,ZN3,XYZN,XYZN1,XYZN2,XYZN3:REAL;
  295. BEGIN
  296. {takes a point x,y,z and rotates the coordinate system its in and gives the new coordinates
  297. for x,y,z in the old coordinate system's perpective.  x1,y1,z1 is new z axis vector, x2,y2,z2 is
  298. new x axis vector, t1,t2,t3 is new coordinates for old systems origin, xm,ym,zm are amplitude
  299. adjusts for x,y,z.  Mode 1 then new point is loaded into r:array[1..3] of real in rectangular
  300. coordinates else r is loaded with spherical coordinates with r[1]=ro,r[2]=theta,r[3]=phi}
  301. IF (X=0)AND(Y=0)AND(Z=0) THEN
  302. BEGIN
  303. U[1]:=X*XM;
  304. U[2]:=Y*YM;
  305. U[3]:=Z*ZM;
  306. GOTO 3;
  307. END;
  308. LEN:=VLENGTH(X,Y,Z,0,0,0);
  309. CROSSPRODUCT(X1,Y1,Z1,X2,Y2,Z2,0,0,0,U3);
  310. YX:=U3[5];
  311. YY:=U3[6];
  312. YZ:=U3[7];
  313. VUNIT(X1,Y1,Z1,0,0,0,TRUE,U4);
  314. ZX:=U4[1];
  315. ZY:=U4[2];
  316. ZZ:=U4[3];
  317. VUNIT(X2,Y2,Z2,0,0,0,TRUE,U4);
  318. XX:=U4[1];
  319. XY:=U4[2];
  320. XZ:=U4[3];
  321.  
  322. H10[1,1]:=XX;
  323. H10[1,2]:=XY;
  324. H10[1,3]:=XZ;
  325. H10[1,4]:=X;
  326. H10[2,1]:=YX;
  327. H10[2,2]:=YY;
  328. H10[2,3]:=YZ;
  329. H10[2,4]:=Y;
  330. H10[3,1]:=ZX;
  331. H10[3,2]:=ZY;
  332. H10[3,3]:=ZZ;
  333. H10[3,4]:=Z;
  334. XN1:=H10[1,4]*((H10[2,2]*H10[3,3])-(H10[2,3]*H10[3,2]));
  335. XN2:=-H10[1,2]*((H10[2,4]*H10[3,3])-(H10[2,3]*H10[3,4]));
  336. XN3:=H10[1,3]*((H10[2,4]*H10[3,2])-(H10[2,2]*H10[3,4]));
  337. XN:=XN1+XN2+XN3;
  338. YN1:=H10[1,1]*((H10[2,4]*H10[3,3])-(H10[2,3]*H10[3,4]));
  339. YN2:=-H10[1,4]*((H10[2,1]*H10[3,3])-(H10[2,3]*H10[3,1]));
  340. YN3:=H10[1,3]*((H10[2,1]*H10[3,4])-(H10[2,4]*H10[3,1]));
  341. YN:=YN1+YN2+YN3;
  342. ZN1:=H10[1,1]*((H10[2,2]*H10[3,4])-(H10[2,4]*H10[3,2]));
  343. ZN2:=-H10[1,2]*((H10[2,1]*H10[3,4])-(H10[2,4]*H10[3,1]));
  344. ZN3:=H10[1,4]*((H10[2,1]*H10[3,2])-(H10[2,2]*H10[3,1]));
  345. ZN:=ZN1+ZN2+ZN3;
  346. XYZN1:=H10[1,1]*((H10[2,2]*H10[3,3])-(H10[2,3]*H10[3,2]));
  347. XYZN2:=-H10[1,2]*((H10[2,1]*H10[3,3])-(H10[2,3]*H10[3,1]));
  348. XYZN3:=H10[1,3]*((H10[2,1]*H10[3,2])-(H10[2,2]*H10[3,1]));
  349. XYZN:=XYZN1+XYZN2+XYZN3;
  350. IF ((ABS(XYZN))<0.000000001) THEN
  351. BEGIN
  352. U[1]:=X*XM;
  353. U[2]:=Y*YM;
  354. U[3]:=Z*ZM;
  355. END
  356. ELSE
  357. BEGIN
  358. U[1]:=(XN/XYZN)*XM;
  359. U[2]:=(YN/XYZN)*YM;
  360. U[3]:=(ZN/XYZN)*ZM;
  361. END;
  362. 3:
  363. IF MODE=1 THEN
  364. BEGIN
  365. U1[1]:=U[1]+T1;
  366. U1[2]:=U[2]+T2;
  367. U1[3]:=U[3]+T3;
  368. END
  369. ELSE
  370. BEGIN
  371. U1[1]:=U[1]+T1;
  372. U1[2]:=U[2]+T2;
  373. U1[3]:=U[3]+T3;
  374. SPHERICAL_XYZ(U1[1],U1[2],U1[3],2,U2);
  375. END;
  376. IF MODE=1 THEN
  377. BEGIN
  378. R1(R)[1]:=U1[1];
  379. R1(R)[2]:=U1[2];
  380. R1(R)[3]:=U1[3];
  381. END
  382. ELSE
  383. BEGIN
  384. R1(R)[1]:=U2[1];
  385. R1(R)[2]:=U2[2];
  386. R1(R)[3]:=U2[3];
  387. END;
  388. END;
  389.  
  390.  
  391. PROCEDURE BAIRSTOW(VAR A,B;EPSILON:REAL;IT1,DEGREE:INTEGER);
  392. LABEL 1,2,3,4,5,6,8;
  393. TYPE
  394. A1=ARRAY[1..200] OF REAL;
  395. B1=ARRAY[1..200,1..2] OF REAL;
  396. VAR
  397. U,U1,V,V1,P,Q,W,R,DENOM,Z,Z1,SUM,RAD:REAL;
  398. F,F1,F2,C1,C2,IT,N:INTEGER;
  399. C:ARRAY[1..200] OF REAL;
  400. D:ARRAY[1..200] OF REAL;
  401. E:ARRAY[1..200] OF REAL;
  402. BEGIN
  403. {bairstow takes a polynomial of to degree 200 and solves for the roots.  if polynomial is
  404. ax^2 + bx + c=0 then first it is normalized to x^2 + (b/a)x + (c/z)=0, then the first
  405. element in a , a[1]=(b/a)  and a[2]=(c/a).  The highest degree coefficient for x^2 is always
  406. understood to be one.  epsilon=as low a number you want for accuracy, it1=number of iterations
  407. before routine calls it quits (for example if the root won't converge), and in this example,
  408. degree is 2.  Bairstow will load B with real and imaginary root values.  b[x,1]=real,b[x,2]=imaginary}
  409. FOR F:=1 TO DEGREE DO
  410. BEGIN
  411. E[F]:=A1(A)[F];
  412. END;
  413. C1:=1;
  414. U1:=0;
  415. V1:=0;
  416. N:=DEGREE;
  417. 1:
  418. IF N<1 THEN
  419. BEGIN
  420. FOR F:=1 TO 200 DO
  421. BEGIN
  422. FOR F1:=1 TO 2 DO
  423. BEGIN
  424. B1(B)[F,F1]:=0;
  425. END;
  426. END;
  427. END
  428. ELSE
  429. BEGIN
  430. IF N=1 THEN
  431. BEGIN
  432. P:=-E[1];
  433. Q:=0;
  434. B1(B)[C1,1]:=P;
  435. B1(B)[C1,2]:=Q;
  436. C1:=C1+1;
  437. END
  438. ELSE
  439. BEGIN
  440. IF N=2 THEN
  441. BEGIN
  442. U:=E[1];
  443. V:=E[2];
  444. 4:
  445. P:=-U/2;
  446. RAD:=(SQR(U))-(4*V);
  447. IF RAD<0 THEN
  448. BEGIN
  449. RAD:=-RAD;
  450. Q:=(SQRT(RAD))/2;
  451. B1(B)[C1,1]:=P;
  452. B1(B)[C1,2]:=Q;
  453. N:=N-1;
  454. C1:=C1+1;
  455. Q:=-Q;
  456. 2:
  457. B1(B)[C1,1]:=P;
  458. B1(B)[C1,2]:=Q;
  459. N:=N-1;
  460. C1:=C1+1;
  461. IF N<=0 THEN GOTO 8;
  462. FOR F2:=1 TO N DO
  463. BEGIN
  464. E[F2]:=C[F2];
  465. END;
  466. GOTO 1;
  467. END
  468. ELSE
  469. BEGIN
  470. Q:=(SQRT(RAD))/2;
  471. W:=P;
  472. R:=Q;
  473. P:=P+Q;
  474. Q:=0;
  475. B1(B)[C1,1]:=P;
  476. B1(B)[C1,2]:=Q;
  477. N:=N-1;
  478. C1:=C1+1;
  479. P:=W-R;
  480. GOTO 2;
  481. END;
  482. END
  483. ELSE
  484. BEGIN
  485. U:=U1;
  486. V:=V1;
  487. IT:=1;
  488. 5:
  489. C[1]:=E[1]-U;
  490. C[2]:=E[2]-(C[1]*U)-V;
  491. FOR F2:=3 TO N DO
  492. BEGIN
  493. C[F2]:=E[F2]-(C[F2-1]*U)-(C[F2-2]*V);
  494. END;
  495. D[1]:=C[1]-U;
  496. D[2]:=C[2]-(D[1]*U)-V;
  497. C2:=N-1;
  498. FOR F2:=3 TO C2 DO
  499. BEGIN
  500. D[F2]:=C[F2]-(D[F2-1]*U)-(D[F2-2]*V);
  501. END;
  502. IF N=3 THEN GOTO 6;
  503. DENOM:=(D[N-1]*D[N-3])-(SQR(D[N-2]));
  504. IF DENOM=0 THEN GOTO 8;
  505. Z:=(C[N]*D[N-3]-C[N-1]*D[N-2])/DENOM;
  506. GOTO 3;
  507. 6:
  508. DENOM:=D[N-1]-(SQR(D[N-2]));
  509. IF DENOM=0 THEN GOTO 8;
  510. Z:=(C[N]-(C[N-1]*D[N-2]))/DENOM;
  511. 3:
  512. Z1:=((D[N-1]*C[N-1])-(D[N-2]*C[N]))/DENOM;
  513. U:=U+Z;
  514. V:=V+Z1;
  515. SUM:=ABS(Z)+ABS(Z1);
  516. IF SUM<EPSILON THEN GOTO 4;
  517. IT:=IT+1;
  518. IF IT>IT1 THEN GOTO 8;
  519. GOTO 5;
  520. END;
  521. END;
  522. END;
  523. 8:
  524. END;
  525.  
  526. PROCEDURE MULTIPLYCR(X,Y,X1,Y1:REAL;VAR Z);
  527. TYPE
  528. Z1=ARRAY[1..2] OF REAL;
  529. VAR
  530. A,B:REAL;
  531. Z2:ARRAY[1..2] OF REAL;
  532. Z3:ARRAY[1..2] OF REAL;
  533. BEGIN
  534. {multiplycr takes two rectangular complex numbers x,y and x1,y1 and stores their product
  535. in z:array[1..2] of real;,z[1]=real,a[2]=imaginary}
  536. COMPLEXCONVERT(X,Y,0,Z2);
  537. COMPLEXCONVERT(X1,Y1,0,Z3);
  538. A:=Z2[1]*Z3[1];
  539. B:=Z2[2]+Z3[2];
  540. COMPLEXCONVERT(A,B,1,Z2);
  541. Z1(Z)[1]:=Z2[1];
  542. Z1(Z)[2]:=Z2[2];
  543. END;
  544.  
  545. PROCEDURE DEVIDECR(N1,N2,D1,D2:REAL;VAR R);
  546. TYPE
  547. R1=ARRAY[1..2] OF REAL;
  548. VAR
  549. A,B:REAL;
  550. R2:ARRAY[1..2] OF REAL;
  551. R3:ARRAY[1..2] OF REAL;
  552. BEGIN
  553. {devidecr takes complex numbers(rect) n1,n2 and d1,d2 and devides n by d storing in r:array[1..2] of real;}
  554. COMPLEXCONVERT(N1,N2,0,R2);
  555. COMPLEXCONVERT(D1,D2,0,R3);
  556. A:=R2[1]/R3[1];
  557. B:=R2[2]-R3[2];
  558. COMPLEXCONVERT(A,B,1,R2);
  559. R1(R)[1]:=R2[1];
  560. R1(R)[2]:=R2[2];
  561. END;
  562.  
  563. PROCEDURE ADDSUBCP(A1,A2,B1,B2:REAL;MODE:INTEGER;VAR R);
  564. TYPE
  565. R1=ARRAY[1..2] OF REAL;
  566. VAR
  567. A,B:REAL;
  568. R2:ARRAY[1..2] OF REAL;
  569. R3:ARRAY[1..2] OF REAL;
  570. BEGIN
  571. {addsubcp takes complex numbers in polar form and adds them for mode=1 or subtracts b from a and
  572. stores in r:array[1..2] of real; with r[1]=magnitude and r[2]=angle}
  573. COMPLEXCONVERT(A1,A2,1,R2);
  574. COMPLEXCONVERT(B1,B2,1,R3);
  575. IF MODE=1 THEN
  576. BEGIN
  577. A:=R2[1]+R3[1];
  578. B:=R2[2]+R3[2];
  579. COMPLEXCONVERT(A,B,0,R2);
  580. R1(R)[1]:=R2[1];
  581. R1(R)[2]:=R2[2];
  582. END
  583. ELSE
  584. BEGIN
  585. A:=R2[1]-R3[1];
  586. B:=R2[2]-R3[2];
  587. COMPLEXCONVERT(A,B,0,R2);
  588. R1(R)[1]:=R2[1];
  589. R1(R)[2]:=R2[2];
  590. END;
  591. END;
  592.  
  593. PROCEDURE POWERCR_CP(A1,A2,B1,B2:REAL;MODE:INTEGER;VAR R);
  594. TYPE
  595. R1=ARRAY[1..2] OF REAL;
  596. VAR
  597. A,B,C,D,E:REAL;
  598. R2:ARRAY[1..2] OF REAL;
  599. R3:ARRAY[1..2] OF REAL;
  600. BEGIN
  601. {if mode is 1 then powercr_cp takes rect complex numbers a and b and raises a to power of b else
  602. powercr_cp takes to polar comples numbers and raises a to b, results stored, as usual, in r}
  603. IF MODE=1 THEN
  604. BEGIN
  605. COMPLEXCONVERT(A1,A2,0,R2);
  606. A:=POWER1(R2[1],B1,0);
  607. B:=EXP((-R2[2])*B2);
  608. C:=LN(R2[1]);
  609. D:=C*B2;
  610. E:=D+R2[2]*B1;
  611. COMPLEXCONVERT(A*B,E,1,R3);
  612. R1(R)[1]:=R3[1];
  613. R1(R)[2]:=R3[2];
  614. END
  615. ELSE
  616. BEGIN
  617. A:=POWER1(A1,B1,0);
  618. B:=EXP((-A2)*B2);
  619. C:=LN(A1);
  620. D:=C*B2;
  621. E:=D+A2*B1;
  622. R1(R)[1]:=A*B;
  623. R1(R)[2]:=E;
  624. END;
  625. END;
  626.  
  627. PROCEDURE RAYS_XYZ(X,Y,Z,X1,Y1,Z1,X2,Y2,Z2,XM1,YM1,ZM1:REAL;COUNT:INTEGER;VAR RIN,ROUT);
  628. TYPE
  629. RIN1=ARRAY[1..100,1..3] OF REAL;
  630. ROUT1=ARRAY[1..100,1..3] OF REAL;
  631. VAR
  632. A:INTEGER;
  633. H2:ARRAY[1..3] OF REAL;
  634. BEGIN
  635. {rays_xyz takes an array of xyz numbers,in rin, up to a 100, and rotates and displaces them
  636. and places the results in rout,  x,y,z is new coordinates for origin, x1,y1,y3 new vector for
  637. z axis, x2,y2,z2 new x axis vector, xm1,ym1,zm1 amplitudes for points, count= number of points to be
  638. processed, rin =incoming array of points, rout=outgoing array of points}
  639. FOR A:=1 TO COUNT DO
  640. BEGIN
  641. XYZROTTRANS(RIN1(RIN)[A,1],RIN1(RIN)[A,2],RIN1(RIN)[A,3],X1,Y1,Z1,X2,Y2,Z2,X,Y,Z,XM1,YM1,ZM1,1,H2);
  642. ROUT1(ROUT)[A,1]:=H2[1];
  643. ROUT1(ROUT)[A,2]:=H2[2];
  644. ROUT1(ROUT)[A,3]:=H2[3];
  645. END;
  646. END;
  647.  
  648. END.
  649.  
  650.  
  651.  
  652.  
  653.  
  654.  
  655.